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We present a study of black hole threshold phenomena for a self-gravitating, massive complex 
scalar field in spherical symmetry. We construct Type I critical solutions dynamically by tuning a 
one-parameter family of initial data consisting of a boson star and a massless real scalar field. The 
massless field is used to perturb the boson star via a purely gravitational interaction which results in 
a significant transfer of energy from the massless field to the massive field. The resulting (unstable) 
critical solutions, which display great similarity with unstable boson stars, persist for a finite time 
before either dispersing most of the mass to infinity (leaving a diffuse remnant) or forming a black 
hole. To further the comparison between our critical solutions and boson stars, we verify and extend 
, the linear stability analysis of Gleiser and Watkins [M. Gleiser and R. Watkins, Nucl. Phys. B319 

■ 733 (1989)] by providing a method for calculating the radial dependence of boson star quasinormal 

modes of nonzero frequency. The frequencies observed in our critical solutions coincide with the 
^SJ , mode frequencies obtained from perturbation theory, as do the radial profiles of many of the modes. 

For critical solutions less than 90% of the maximum boson star mass M max ~ 0.633M£,/to, the 
existence of a small halo of matter in the tail of the solution distorts the profiles which otherwise 
agree very well with unstable boson stars. These halos appear to be artifacts of the collision between 
the original boson star and the massless field, and do not appear to belong to the true critical 
solutions, which are interior to the halos and which do in fact correspond to unstable boson stars. 
It appears that unstable boson stars are unstable to dispersal ("explosion") in addition to black 
hole formation, and given the similarities in macroscopic stability between boson stars and neutron 
stars, we suggest that those neutron star configurations at or beyond the point of instability may 
likewise be unstable to explosion. 
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I. INTRODUCTION 



Over the past decade, detailed studies of models of gravitational collapse have revealed that the threshold of black 
5jT)i hole formation is generically characterized by special, "critical" solutions. The features of these solutions are known 
as "critical phenomena," and arise in even the simplest collapse models, such as a model consisting of a single, real, 
massless scalar field, minimally coupled to the general relativistic field in spherical symmetry [jl]. Although we present 
a brief overview of black hole critical phenomena here, we suggest that interested readers consult the excellent reviews 
by Gundlach |^,|| for many additional details. 

Black hole critical solutions can be constructed dynamically via simulation, i.e. via solution of the full time- 
dependent PDEs describing the particular model, by considering one-parameter families of initial data which are 
required to have the following "interpolating" property: for sufficiently large values of the family parameter, p, the 
evolved data describes a spacetime containing a black hole, whereas for sufficiently small values of p, the matter-energy 
in the spacetime disperses to large radii at late times, and no black hole forms. For any such family, there will exist a 
critical parameter value, p = p*, which demarks the onset, or threshold, of black hole formation. To date at least, it 
has invariably turned out that the solutions which appear in the strongly-coupled regime of the calculations (i.e. the 
critical solution), are almost totally independent of the specifics of the particular family used as a generator. In fact, 
the only initial-data dependence which has been observed so far in critical collapse occurs in models for which there 
is more than one distinct black-hole-threshold solution. In this sense then, black hole critical solutions are akin to, 
for example, the Schwarzschild solution, which can be formed through the collapse of virtually any type and/or shape 
of spherically distributed matter. In particular, like the Schwarzschild solution, black hole critical solutions possess 
additional symmetry (beyond spherical symmetry) which, to date, has either been a time-translation symmetry, in 
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which the critical solution is static or periodic, or a scale-translation symmetry (hometheticity) , in which the critical 
solution is either continuously or discretely self-similar (CSS or DSS). 

However, in clear contrast to the Schwarzschild solution, black hole threshold solutions are, by construction, unstable. 
Indeed, after seminal work by Evans and Coleman [Q and by Koike et al || , we have come to understand that critical 
solutions are in some sense minimally unstable, in that they tend to have precisely one unstable mode in linear 
perturbation theory. Thus letting p — > p* amounts to minimizing or "tuning away" the initial amplitude of the 
unstable mode present in the system. 

As already suggested, two principal types of critical behavior have been seen in black hole threshold studies; which 
type is observed depends, in general, upon the type of matter model and the initial data used — as mentioned, some 
models exhibit both types of behavior. Historically, one of us termed these Type I and Type II solutions, in a 
loose analogy to phase transitions in statistical mechanics, but at least at this juncture, we could equally well label 
the critical solutions by their symmetries (i.e. static/periodic or CSS/DSS) . For Type I solutions, there is a finite 
minimum black hole mass which can be formed, and, in accord with their static/periodic nature, there is a scaling 

law, r ~ —7 In \p—p* I , relating the lifetime, r, of a near-critical solution to the proximity of the solution to the critical 
point. Here 7 is a model-specific exponent which is the reciprocal of the real part of the eigenvalue associated with 
the unstable mode. On the other hand, Type II critical behavior — less relevant to the current study — is characterized 
by arbitrarily small black hole mass at threshold, and critical solutions which are generically self-similar. 

The direct construction, or simulation, of critical solutions, has thus far been performed almost exclusively within 
the ansatz of spherical symmetry. In the spherical case one must couple to at least one matter field for non-trivial 
dynamics, and spherically symmetric critical solutions for a considerable variety of models have now been constructed 
and analyzed. In addition to the massless scalar case mentioned above, these include solutions containing a perfect 
fluid J0,|| , a scalar non- Abelian gauge field J 7] , and particularly germane to the current work, a massive real scalar 
field H. The work of Abrahams and Evans which considered axisymmetric critical collapse of gravitational waves, 
remains notable for being the only instance of a reasonably well-resolved non-spherical critical solution p0{ . 

Our current interest is a critical-phenomena- inspired study of the dynamics associated with "boson stars" [PHhSl , 
a class of equilibrium solutions to the Einstein-Klein Gordon system for massive complex fields, which are supported 
against gravitational collapse by the effective pressure due to the dispersive nature of a massive Klein-Gordon field. 
(For extensive reviews on the subject of boson stars, see Jetzer [[bi} or Mielke and Schunck fig] .) We know from the 
studies by Gleiser and Watkins and by Lee and Pang jlTj , that there exists a critical value of the central density 
which marks the transition between boson stars which are stable with respect to infinitesimal radial perturbations, 
and those which are unstable. The dynamical simulations of Seidel and Suen |l8| revealed scenarios in which a boson 
star on the unstable branch would either form a black hole or radiate scalar material and form a boson star on the 
stable branch. Their study is extended in this paper, in which we consider dynamical changes to the geometry of a 
boson star which are large enough to bring it to the threshold of black hole formation. 

As already mentioned, another paper closely related to this work is that of Brady et al. which described a 
dynamical study of critical solutions of a massive real scalar field. Those authors demonstrated scenarios in which 
black holes could be formed with arbitrarily small mass (Type II transitions), and those in which the black holes formed 
had a finite minimum mass (Type I transitions) . The boundary between these regimes seemed to be the relative length 
scale of the pulse of initial data compared to the Compton wavelength associated with the boson mass. Initial data 
which was "kinetic energy dominated" evolved in a manner essentially similar to the evolution of a massless scalar 
field. Initial data pulses having widths larger than the length scale set by the boson mass were "potential dominated," 
providing a characteristic scale for the formation of the critical solutions. Brady et al. found that the resulting Type 
I critical solutions corresponded to a class of equilibrium solutions discovered by Seidel and Suen which are 
called "oscillating soliton stars." These soliton stars share many characteristic with the complex-valued boson stars, 
such as the relationship between the radius and mass of the star. Both types of "stars" have a maximum mass, 
and show the same overall behavior as "real" (fermion) stars in terms of the turn-over in their respective stability 
curves. Interestingly, although the soliton stars are not static — they are periodic (or quasi-periodic) — many of the 
same macroscopic properties seen in fluid stars are still observed. 

In this paper, we construct critical solutions of the Einstein equations coupled to a massive, complex scalar field 
dynamically, by simulating the implosion of a spherical shell of massless real scalar field around an "enclosed" boson 
star. The massless field implodes toward the boson star and the two fields undergo a (purely gravitational) "collision." 
The massless pulse then passes through the origin, explodes and continues to r — > 00, while the massive complex 
(boson star) field is compressed into a state which ultimately either forms a black hole or disperses. We can thus 
play the "interpolation game" using initial data which result in black hole formation, and initial data which give 
rise to dispersal: specifically, we vary the initial amplitude of the massless pulse to tune to a critical solution. We 
analyze the black hole threshold solutions obtained in this manner, and discuss the similarities between our critical 
solutions for the self-gravitating complex massive scalar field and boson stars on the unstable branch. To further this 
discussion, we extend the work of Gleiser and Watkins um and compare the results of the simulations with those of 
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linear perturbation theory. 

The layout of the remainder this paper is as follows: In Section II, we describe the mathematical basis for our 
numerical simulations. In Section III, we present results from our simulations, in which the Type I character of 
the critical solutions is demonstrated, along with the close similarities one finds between the features of the critical 
solutions and those of boson stars. In most of the critical solutions we find a halo of mass near the outer edge of 
the solution which is not a feature of boson star equilibrium data. Inside this halo, however, the critical solutions 
match the boson star profiles very well. In Section IV, we give a synopsis of our linear stability analysis of boson star 
quasinormal modes, from which we obtain the boson star mode frequencies as functions of the central value of the 
modulus of the complex field. Section V concerns the radial profiles of the perturbative modes per se, and includes 
a comparison of the mode shapes and frequencies obtained from perturbation theory with our simulation data. The 
modes obtained by these two different methods agree well with each other, although the additional oscillatory mode 
in our simulation data is only shown to agree with the corresponding boson star mode in terms of the oscillations in 
the metric and not in the field (possibly as an artifact of our simplistic approach to extracting this mode from the 
simulation). In Section V we provide further discussion regarding the properties of the halos surrounding the critical 
solutions. 

Conclusions in Section VI are followed by appendices giving tables of mode frequencies versus the central field value 
of the boson star, details about our finite difference code, and details of our linear stability analysis, which includes 
a description of our algorithm for finding the frequencies of boson star modes. 



II. SCALAR FIELD MODEL 



A boson star is described by a complex massive scalar field 0, minimally coupled to gravity as given by general 
relativity. We work solely within the context of classical field theory, and choose units in which G and c are unity. 
Furthermore, since all lengths in the problem can be scaled by the boson mass m ]l3fl , we choose m — 1. To the usual 
boson star model, we add an additional, massless real scalar field, 03, which is also minimally coupled to gravity. This 
additional scalar field will be used to dynamically "perturb" the boson star. 

The equations of motion for the system are then the Einstein equation and Klein-Gordon equations: 

Gab = Rab - \g ab R = 8^ (Tf b (0) + T* b (fo)) (2.1) 

□ 0-m 2 = O (2.2) 

□ 3 = (2.3) 

where 

87tT q c & (0) = d a <f>db^* + d Q 0*d b - gab (d c 0d c 0* + m 2 |0| 2 ) , (2.4) 



87^(03) = 2d a 3 d fc 03 - 5afcd c 3 d c 03, (2.5) 

and □ is the DAlembertian operator. While more general potentials in ( |2.2| ) have been employed recently pO| , pT| , 
we will restrict our discussion to the simplest case, i.e. merely the m 2 2 potential. We also stress that the complex 
scalar field, 0, and the massless, real scalar field, 3 are coupled only through gravity — in particular we do not include 
any interaction potential V/(0, 3 ). 

Restricting our attention to spherical symmetry, we write the most general spherically-symmetric metric using 
Schwarzschild-like "polar-areal" coordinates 

ds 2 = -a 2 (t, r)dt 2 + a 2 (t, r)dr 2 + r 2 dil 2 , (2.6) 

and generally make use of the "3+1" formalism of Arnowitt, Deser and Misner |2^] which regards spacetime as a 
foliation of spacelike hypersurfaces parameterized by t. 

We write the (spherically-symmetric) complex field, 0(i,r), in terms of its components 

0(t,r) =0 1 (*,r)+i0 2 (t,r) (2.7) 



3 



where (f>x(t, r) and 4>2(t,r) are each real. Again, since our model includes no self-interaction (anharmonic) potential 
for the complex field, </>i and <f>2 are only coupled through the gravitational field. 
We then define 



$i(t,r) = # 



(2.8) 



Hi (t,r) = -0i 
a 



U 2 {t,r) = -(f>2, 
a 



(2.9) 



n 3 (*,r) = -0 3 . 
a 



where ' = <9/dr and ' = 9/<?i. 

With these definitions, the equations we solve are the Hamiltonian constraint, 



a 
a 



1-a 2 
2r 



- [il 2 + n 2 2 + n 3 2 + $ x 2 + $ 2 2 + $ 3 2 + a > 



(where 111 2 should be read as (111) 2 ), the slicing condition 

a' _ a 2 - f a' 
a 

and the Klein- Gordon equations, 



+ --a 2 r{<t>t 2 + <p 2 % 
r a 



( T 2 Ot 



(2.10) 



(2.11) 



(2.12) 



(2.13) 



where k = 1, 2, 3 and <!)3fc is a Kronecker delta used to denote the fact that 03 is a massless field. 

We also have equations which are used to update the spa tial gr adie nts of the scalar fields, as well as the scalar fields 
themselves. These follow directly from the definitions (2.S) and (|2.9|): 



(2.14) 



$ k dr 



(2.15) 



Equations (pll])-(|2.15|) are solved numerically using the second order finite difference method described in Appendix 
B. 

Initial conditions for our simulations are set up as follows. First, initial data for the massive field are constructed 
from the boson star ansatz 



0(t,r) = (j) {r)e~ 



(2.16) 



where we let <po{r) be real. Substitution of this ansatz into the full set of equations ( 2.11 )-( 2.15 ), yields a system 
of ordinary differential equations (ODEs), whose solution, for a given value of the central field modulus, is found by 
"shooting," as described in frl^] . Once the boson star data is in hand, we add the perturbing massless field by freely 
specifying $3 and II3. At this point, all matter quanti ties h ave b een s pecified ; the initial geometry, et(0, r) and a(0, r) 
is then fixed by the constraint and slicing equations ( 2.11 ) and ( 2.12 ). 

In relating the simulation results which follow, it is useful to consider the individual contributions of the complex 
and real fields to the total mass distribution of the space-time, in order that we can meaningfully and unambiguously 
discuss, for example, the exchange of mass-energy from the real, massless field to the massive, complex field. By 
Birchoff's theorem, in any vacuum region, the mass enclosed by a sphere of radius r at a given time t is given by the 
Schwarzschild-like mass aspect function M(t,r) — r(l — l/a 2 )/2. However, at locations occupied by matter, M(t,r) 
cannot necessarily be usefully interpreted as a "physical" mass. In polar-areal coordinates, the mass aspect function 
is related to the local energy density p(t, r) by 



dM(t, r) 
dr 



(2.17) 
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with p(t, r) given in our case by 



p{t, r) = —j [n! 2 + n 2 2 + $! 2 + $ 2 2 + a 2 (0! 2 + 2 2 )] + — [n 3 2 + +$ 3 2 ] . (2.18) 



Here, we have explicitly separated the contributions from the complex and real fields. Since dM/dr is given by a 
linear combination of the contributions from each field, we can decompose dM/dr so that, in instances where there is 
no overlap in the supports of the distinct fields, we can unambiguously refer to the mass due to one field or the other. 
That is, we can refer to the individual contributions of each field to the total mass as being physically meaningful 
masses in their own rights. Then, by integrating the contribution of each field to dM/dr over some range of radius 
(fmin ■ • •''max), (where there is some region of vacuum starting at r m i n and extending inward, and some region of 
vacuum starting at r >— r max and extending outward), and demanding that none of the other type of field is present 
in the domain of integration, we can obtain a measure of the mass due to each field. 

Motivated by such considerations, we define an energy density for the complex field, pc, as 

Pc(t, r) = [ ni2 + ^ 2 + + $2 ' + fl2 + ^ 2 )] ' ( 2 - 19 ) 
with a corresponding mass aspect function, Mc{t,r), given by 

M c {t,r)= ( r 2 p c dP. (2.20) 



Similarly, the energy density due to the real field is defined as 



2a 2 

with a corresponding mass aspect function, Mn(t,r) given by 



p R (t,r)^ 7 ^[U 3 2 + +^ 3 2 ], (2.21) 



M R (t,r) = / r*p R df. 
Jo 

We again emphasize that in regions where the supports of the different fields overlap (and in non-vacuum regions 
in general) it may not be possible to ascribe physical meaning to the individual mass aspect functions defined above. 
(However, even in such instances, these functions are still useful diagnostics.) Most importantly, where the supports 
of the fields do overlap, and only in these regions, it is possible for mass-energy to be exchanged from one scalar field 
to the other — through the gravitational field — while the sum Mc + Mr — M (measured in an exterior vacuum region) 
is conserved. The quantities given above allow us to measure this exchange of mass by looking at the profiles Mc(t, r) 
and Mfl(t, r) before and after a time when the fields are interacting. This is shown in the next section. 

As a further consideration, we point out that the U(l) symmetry of the complex field implies that there is a 
conserved Noether current, J 11 , given by 

^ = ^g^m,? - <f>*d v <t>). (2.22) 

071" 

The corresponding conserved charge or "particle number" N is 

N -- 

We may also wish to regard N as a function of t and r by integrating the above function from zero to some finite 
radius, in which case 

2*^ = ^(^-11^!). (2.23) 

Some authors have considered the difference Mc — mN to be a sort of "binding energy" of the complex field juj , 
however this quantity does not correspond to any transition in the stability of boson stars, and we have not found it 
to be very useful in understanding the dynamics 

of our simulations. 

Finally, following Seidel and Suen Jlq] , we define a radius Rg5(t,r) for the boson star implicitly by Mc\r 95 — 
0.95 Mc\ r — too ■ Alternatively, we will also consider a radius Res(t,r) which encloses (1 — e _1 ) ~ 63% of Mc 
and which will include the "bulk" of a boson star but will neglect the "tail" . 
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III. SIMULATION RESULTS 



We choose the initial data for the complex field to be a boson star at the origin, with a given central density 4>o(Q). 
For the massless field </>3(0,r), we choose one of the families in Table |l| We generate critical solutions by tuning the 
amplitude A of 03(0, r) (holding the position ro and width A constant) using a bisection search, until the resulting 
solution is arbitrarily close (i.e. within some specified precision) to the point of unstable equilibrium between dispersal 
and black hole formation. 

Figure [l] shows a series of snapshots from a typical simulation in which the parameter p (p = A) , is slightly below 
the critical value p*, for a boson star on the stable branch with a mass of M ~ 0.59AfJ, ; /m (where Mpi is the Planck 
mass). The shell of massless field, a member of initial data Family I, implodes through the boson star and explodes 
back out from the origin, and the gravitational interaction between the fields forces the boson star into a new state, a 
"critical solution." We see from this animation, and from Figure 3, that dispersal from the critical state does not mean 
that the boson star returns to its original stable configuration, but rather that the star becomes strongly disrupted 
and "explodes." That is to say, if we were to follow the evolution beyond t = 475, the massive field would continue 
to spread toward spatial infinity. At some late time, after a large amount of scalar radiation has been emitted, the 
end state would probably be a stable boson star with very low mass. 



TABLE I. Families of initial data. For both fa milie s, the i nitial data, <f>(0, r) = </>i(0,r) + zcfe(0, r~), f or the massive complex 
field is given by a boson star, obtained by solving (2.11 )— ( 2.13 ) numerically according to the ansatz (2.16) (with the parameter ui 
found via "shooting"). The initial real massless field profile, </>3(0, r), is given in closed form by the "gaussian" and "kink" initial 
data. For each family, we also choose d t (j>3(0,r) such that the pulse is initially in-going, i.e. 113(0, r) = $3(0, r) + cj>z(0,r)/r. 



Family 



I 

II 



Complex Field <j>\ + i<\>2 
Name Parameters Profile 



Boson Star </>o(0) <t>a{ r ) 
Boson Star 0o(O) <t>o(r) 



Real Field </>3 
Name Parameters Profile 



Gaussian A, tq , A A exp 



.4 



Kink A, r ,A — 1 + tanh . 
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FIG. 1. Evolution of a perturbed boson star with ^o(O) = 0.04 x and mass Mc = 0.59M|> ; /m. This shows contributions 
to dM/dr due to the massive field (solid line) and massless field (dashed line). We start with a stable boson star centered at 
the origin, and a pulse of massless field given by Family I with ro = 30 and A = 8. (We see two peaks in dM/dr of massless 
field because it is only the gradients of cf>3, not (j>3 itself, which contribute to Mn(r,t) for a massless field.) In the evolution 
shown above, the pulse of massless field enters the region containing the bulk of the boson star (t ~ 15), implodes through the 
origin (t ~ 30) and leaves the region of the boson star (t ~ 50). Shortly after the massless pulse passes through the origin, the 
boson star collapses into a more compact configuration, about which it oscillates for a long time before either forming a black 
hole or dispersing. (The case of dispersal is shown here.) Note that the perturbing field <f>z passes through the boson star and 
exits the region containing most of the star, even before the massive field reaches its denser, critical state. Thus the massless 
field is not part of the critical solution per se. Black hole formation (always with a finite black hole ADM mass in our study) 
can take place at times long after the massless pulse has left the neighborhood of the boson star. 



The gravitational interaction between the two fields results in an exchange of energy from the massless field to the 
massive field, as shown in Figure ^. Figure ^ shows some timelike slices through the simulation data, giving a plot of 
the maximum value of a, the value of \<j>\ at the origin, and the radius -R95 as functions of time. These are given to 
help elucidate the point that the critical solution oscillates about some local equilibrium, before dispersing or forming 
a black hole. The lifetime of the critical solution increases monotonically as p — > p*. Figure |] shows that the scaling 
law expected for Type I transitions is exhibited by these solutions. 
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00 20 40 60 80 100 



t. 

FIG. 2. Exchange of energy between the real and complex scalar fields. For this simulation, initial data from Family I was 
used, with 4>o(0) = 0.04 x \/47r, ro = 40 and A = 8. The solid line shows the mass of the complex field, shifted upward on the 
graph by 0.21Af|> ; /m. The long-dashed line shows the mass of the real field, shifted upward by 0.55M Pl /m. The mass AM 
exchanged from the massless field to the massive field in this simulation is nearly 0.0053, or about 2.5% of the mass of the real 
field (9% of the boson star mass). The amount (and percentage) of mass transfer goes to zero as we consider boson star initial 
data approaching the transition to instability (see, e.g. Figure 7). The dotted line near the top of the graph shows the total 
mass enclosed within r = 100. Throughout the simulation, both the total mass M = Mc + Mr and the particle number N (of 
the complex field) are conserved to within a few hundredths of a percent. 




100 200 300 400 500 



FIG. 3. Quantities describing a near-critical solution. Here we show timelike slices through the data shown in Figure |l| an 
evolution that ends in dispersal. Top: Maximum value of the metric function a (for each spacelike hypersurface parameterized 
by t). The local maximum at t ~ 40 is due to the presence of the pulse of massless field. Middle: Central value \4>(t, 0)| of the 
massive field. Bottom: Radius 7?gs which contains 95% of the mass-energy in the complex field. Again, we see evidence that 
after the remaining in critical regime for a while, the star can "explode," leaving a diffuse remnant with low mass. 
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FIG. 4. Lifetime r of a typical set of near-critical solutions vs. In \p — p*\. We use initial data from Family I. The lifetime 
of the critical solution obeys a simple scaling relation. Using super-critical solutions, we measure r to be the time from t — 
until black hole formation occurs. The relationship shown in the graph can be described by r = — 7ln|p — p*\, where for the 
data shown in this graph, 7 ~ 9.2 The value of 7 can be related to the imaginary part of the Lyapunov exponent a of the 
unstable mode (~ e i<Tt ) by Qf(o-) = I/7 ~ 0.109. This value is the same as that obtained from a linear perturbation analysis of 
the specific boson star to which we believe this configuration is asymptoting (See Section V). We note that in the limit p — > p* , 
the mass of the black hole formed is finite (and close to the mass of the progenitive unstable boson star) , i.e. the system 
exhibits Type I critical behavior. 
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FIG. 5. Mass vs. radius for equilibrium configurations of boson stars (solid line), initial data for the complex field (triangles), 
and critical solutions (squares) . Arrows are given to help match initial data with the resulting critical solutions. Points on the 
solid line to the left of the maximum mass M max ~ 0.633Mj> ; /m correspond to unstable boson stars, whereas those to the right 
of the maximum correspond to stable stars. If one takes time averages of properties such as mass, central density [0(t,O)| and 
radius R95 during the critical regime, one finds values which match the profile of a boson star on the unstable branch. The 
squares show the time average of each critical solution during the oscillatory phase. Graph (a) shows mass M versus Tigs the 
radius containing 95% of M, whereas graph (b) shows M versus the radius containing (1 — e _1 ) M. The agreement between 
the critical solutions and boson stars shown in graph (a) deteriorates with decreasing mass, however the comparison shown in 
graph (b), which neglects the "tail" of the critical solutions and boson stars, shows much better agreement for all masses. (We 
show the tail region in Figure 6|.) In this simulation the massive field radiates only a small amount due to the perturbation 
by the massless field, and so the stable boson star is essentially driven to "pop" across the stability curve by the impinging 
massless pulse. 
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Figure |5| shows the mass vs. radius for some critical solutions along with the equilibrium curve for boson stars. We 
notice that there are great similarities, at least for relatively high mass configurations, between the critical solutions 
and unstable boson stars in the ground state. (We do not perform studies involving boson stars with much lower 
masses, because of the dynamic range required for the spatial resolution of the finite difference code. Also, for a given 
numerical error tolerance, such low-mass critical solutions have much shorter lifetimes than larger-mass solutions and 
thus it is more difficult to obtain average properties of such short-lived critical solutions.) When we include nearly all 
of the complex-scalar mass in our comparisons, as shown in Figure |^(a), we see that the time-averaged properties of 
the critical solutions with lower masses, i.e. those further from the transition to instability, deviate from the curve 
of equilibrium configurations, and that the deviation increases as mass decreases. When we consider only the bulk of 
the boson star, however, we see very good agreement between the dynamically generated critical solutions and the 
unstable boson stars, computed from the static ansatz, as shown in Figure ^|(b). The comparison between low-mass 
critical solutions and boson stars, shown in Figure |[ can be further illuminated by looking at a profile of the mass 
distribution as shown in Figure M. 
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FIG. 6. Comparison of highly unstable (low-mass) critical solution and boson star. Squares show a critical solution resulting 
from a boson star having (j)o(0) = 0.26 x \fin. (The data has been reduced for graphing purposes; the actual spatial resolution 
in the simulation is four times finer than that shown in the figure.) The solid line shows a "best fit" (unstable) boson star we 
constructed by finding the time average of \<p(t, 0)| in the critical solution and using this as the value for <j>o(0) in the ODE 
integration routine which solves for the equilibrium (boson star) solutions. We see that there is a small halo near the outer edge 
of the solution (r = 8) . The halo has the same relative magnitude when viewed in terms of the particle number distribution 
dN/dr. We discuss the halo phenomena further in Section V. 



We see that there is a small halo near the outer edge of the solution (r = 8), and it is this which throws off our 
measurement of -R95 used for Figure ||. Despite the effect this has on the measurement of the radius R 95 of the star, 
we can still obtain a good fit of a boson star to the interior of the critical solution in the low-mass regime. We provide 
further discussion of these halos in Section V. in the critical solutions with higher total mass. 

It is also worth noting that the critical solution best corresponds to a boson star in the "ground state," i.e., without 
any nodes in the distribution of the fields <j>\ or (f>2- Boson stars in excited states (i.e., having nodes in </>i and fa) 
have mass distributions which differ significantly from the critical solutions we obtain [20| ] . 

We wish to explain these simulation results in terms the quasi-normal modes of boson stars. Previous work in 
critical phenomena jfj] - , |2^] leads us to surmise that there is a single unstable mode present in the system which 
is excited when the boson star moves into the critical regime. The oscillatory behavior during the critical regime may 
be explainable in terms of the superposition of a stable oscillatory mode with the unstable mode. In the next section, 
we attempt to confirm these hypotheses by means of perturbation theory. 
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IV. BOSON STAR STABILITY STUDY VIA LINEAR PERTURBATION THEORY 



We follow the work of Gleiser and Watkins fig] . For the perturbation calculations, we find it helpful to define the 
following metric functions: 



e v ^ = a 2 
e x ^ = a 2 . 



and to rewrite the complex field 4>{t 1 r) as 

4>(t,r) = [^i(i,r)+z^ 2 (i,r)]e" Mt , 



(4.1) 



where ipi and ip2 are real. (Note that this is a different decomposition of the field 4> than (2/7), the one used in the 
previous sections.) 

In these variables, the equilibrium quantities are 



X(t,r) = X (r) 
v(t,r) = v (r) 
ipi(t,r) = (f> {r) 



(4.2) 
(4.3) 
(4.4) 
(4.5) 



For the perturbation, we expand about the equilibrium quantities by first introducing four perturbation fields- 
SX(t, r), Sv(t,r), 8ipi(t,r) and Sip2(t,r) — and then setting: 



X(t,r) = X (r) + 6X(t,r) 
v(t,r) = v (r) + Sv(t,r) 
^ 1 (t,r) = 0o(r)(l + ^i(t,r)) 
ip2(t,r) = 4>o(r)Sip2(t,r). 



(4.6) 
(4.7) 
(4.8) 
(4.9) 



These expressions are substituted into the relevant evolution and constraint equations (details in Appendix C) , after 
which the resulting system can be reduced to the following system of two coupled second-order ordinary differential 
equations for 5<fii and 8X: 
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To perform the stability analysis (normal- mode analysis), we assume a harmonic time dependence, i.e., 



(4-11) 



Sipi (t,r) = 5tpi(r) e 4 ' 
SX(t,r) = SX(r) e lat . 



Note that ( 4.1C ) and ( 4.11 ) contain only second derivatives with respect to time, and because there are good reasons 
to assume a 2 is purely real p^Jlr| , we only need to determine whether a 2 is positive or negative to determine stability 
or instability, respectively. 

Using the method described in Appendix C, we find the distribution for the squared frequency <Tq of the fundamental 
mode, with respect to 4>q, which is shown in Figure 0. 
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Superposed with the fundamental mode, we may have other modes at higher frequencies. Figure |8| shows the 
relation between first harmonic frequencies and (j>o(0). 
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FIG. 7. Mode frequencies of boson stars: fundamental mode This plot shows a graph of ctqi the squared frequency of the 
fundamental mode, versus the value of 4>q at the origin. Note that, as the inset shows, oq crosses zero when ^>o (0) ~ 0.27, 
which corresponds to a boson star with the maximum possible mass. (The circles show actual values obtained, and the solid 
line simply connects these points.) 
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FIG. 8. Mode frequencies of boson stars: first harmonic mode This plot shows a graph of a\, the squared frequency of the 
first harmonic mode, versus the value of <j>o at the origin. Note that, as the inset shows, <j\ crosses zero when </>o(0) — 1.15, 
which corresponds to the first local minimum on the unstable branch of the mass vs. radius curve (see Figure 4). (The circles 
show actual values obtained, and the solid line simply connects these points.) 

V. COMPARISON OF PERTURBATION ANALYSIS AND SIMULATION DATA 

We wish to compare the results of our perturbation theory calculation with the oscillations of stable boson stars. 
Two differences exist between the conventions used in the perturbation theory calculation and those used in the boson 
star simulation data. The first difference is in the choice of the time coordinate. In the perturbation theory code, we 
choose a lapse of unity at the origin, whereas in the simulations we set the lapse to unity at spatial infinity. Thus we 
have the following mapping from the perturbation theory calculations to the simulations: 
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Pcrturbative 



Simulation 



The other significant difference is in the way the complex field <f>(t, r) is decomposed into constituent real fields. 
Thus we cannot directly compare 4>\ and ip\, for example. We can, however, compare the modulus \<p\ of the field. 
For the simulation data, the perturbation in \<j>\ can be taken directly from (<pf + (j} 2 ,) 1 ^ 2 ■ For the data obtained from 
perturbation theory, the perturbation in \<p\ will be, to first order, (froSipi. 

Before proceeding to the comparisons per se, we wish to point out that determining the unstable mode via numerical 
simulation of the full nonlinear system was very easy to do in comparison to the linear perturbation theory calculations. 



A. Unstable modes 



To measure the unstable mode, we again perform a series of simulations in which we allow a gaussian pulse from 
an addition real, massless Klein-Gordon field to impinge on a stable boson star. 

By tuning the amplitude of this pulse (holding constant the width of the pulse and its initial distance from the 
boson star), we can generate a family of slightly different near-critical solutions depending on the amplitude of the 
initial gaussian pulse, and can tune down the initial magnitude of the unstable mode. By subtracting these slightly 
different near-critical solutions, we obtain a direct measurement of the unstable mode. 

Considering a specific example, we start with a stable boson star which has an initial field value at the origin of 
4>0 (0) = 0.04 x \/47r. By driving it with a gaussian pulse tuned to machine precision, we can cause this stable star 
to become a critical solution which persists for very long times, oscillating about a local equilibrium. The average 
value of |</>(i,0)| is (|</>(i,0)|) ~ 0.463. We measure the unstable mode by subtracting data of a run which contained 
a gaussian pulse with an amplitude that differed by 10 -14 from that of the pulse tuned to machine precision. We 
can then measure the growth factor of the unstable mode by taking the L2 norm of this difference at various times, 
taking the logarithm, and fitting a straight line to it. From this, we obtain a ~ 0.109 i, or a 2 ~ —0.0118. Because of 
the differences in time coordinate between the simulations and perturbation theory calculations, we need to compute 
a 2 /a 2 in order to compare with the perturbation calculations. We find the average value of l/a(i,0) 2 for the times 
listed above to be (l/a(t, 0) 2 ) ~ 3.80, and thus we find a 2 /a 2 ~ —0.0450. We choose to compare these perturbation 
theory results with data from a time in the simulation for which the difference in field values (for the two evolutions 
tuned slightly differently) is A|0(t, 0)| ~ 8.4 x 10~ 13 . We use this value in the perturbation theory solver and arrive 
at a 2 ~ —0.045, in good agreement with the value from the simulation. In Figures p| and [Tc|, we compare the graphs 
of the solutions for the unstable mode. In Figure [0] we show a comparison of the squared freqency values obtained 
from the linear perturbative analysis and those as measured in our simulations. 
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FIG. 9. Fundamental mode of unstable boson star, (a) The solid line shows <f>ofii>i from the perturbation theory calculations. 
The squares shows the difference between \<j>\ for two simulations for which the critical parameter p differs by 10 -14 . (The data 
has been reduced for graphing purposes; the actual spatial resolution in the simulation is four times finer than what is shown 
in the figure.) Differences between the simulation data and perturbation theory results are below 1.1 x 10~ 15 . If a line were 
drawn connecting the squares, it would be indistinguishable, to the eye, from the perturbation theory line. Thus we provide a 
second graph (b) showing the difference of these results, where the scale is relative to the maximum value of 6\tj>\. 
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FIG. 10. Fundamental mode of unstable boson star, (a) The solid line shows the perturbation to the metric function a, 
as found from the perturbation theory calculations. The squares shows the difference between the metric function a for two 
simulations for which the critical parameter p differs by 10 -14 . (In the simulations, the spatial resolution was four times that 
shown in the figure.) (b) A plot of the difference between the mode obtained from the simulation and the mode obtained via 
perturbation theory, where the scale is relative to the maximum value of 8a. 
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FIG. 11. Comparison of squared frequencies/ Lyapunov exponents for unstable modes. The circles show a subset of the 
perturbation theory data as displayed in Figure [?]. The squares show the measurements obtained from our simulations. (The 
solid line simply connects the circles.) We note that the agreement between the two sets is good even for the more unstable, 
low-mass solutions. We also point out that the measurements of our simulations were performed along r = 0, i.e.,, in the 
interior of the halo found in the low-mass solutions, which seems to strengthen the remarks at the end of section III, namely 
that, aside from the halo at the exterior of the critical solution, the critical solutions (of all masses) seem to correspond to 
unstable boson stars. 



B. Oscillatory modes 

We can also look at the oscillatory mode during the critical regime. We study the behavior of such a mode using 
the same technique we used to examine the fundamental mode of the unstable boson star: we subtract the data at 
one instant of time from the data at all other instants. Again, as a specific example, we use the same initial boson 
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star as that used in the previous section. During the critical portion of the evolution, we notice an oscillation period 
of about T ~ 38.4, and thus we obtain a — 2-k/T ~ 0.0261. During this period, the average value of l/a 2 (t, 0) is 
about 3.80, and thus we find a 2 /a 2 ~ 0.102. We take data from a moment in the middle of the oscillation period, and 
subtract it from the data at other times. We can then compare the perturbation theory results with simulation data 
at a local peak of the oscillation. For the local peak we chose at time t — t p , the difference in modulus of the field was 
A\cf)(tp, 0)| ~ 0.0197. Inserting this value into the perturbation theory code, we find a 2 ~ 0.102 for this configuration. 
Thus we again find excellent agreement between the squared oscillation frequencies computed in perturbation theory 
and via simulation. 

In Figures [l^ and [l3| we compare the functions obtained from the perturbation theory calculation with those from 
the simulation. We note that the agreement for the metric functions is very good for all radii, but the agreement in 
the fields begins to decline beyond r = 5. Why do the graphs of \<j)\ not agree well for the first harmonic? This could 
be a consequence of our simplistic method of extracting this mode. While our method of simply subtracting different 
frames has worked well for our test cases of oscillations of stable boson stars, the first harmonic of the unstable star 
has a higher frequency and thus our graph could be subject to sampling effects. A better method would be to perform 
a Fourier transform in time for each grid point, and construct the high er harmonics in the field accordingly. There 



may be a simple resolution to the discrepancy in the graphs of \<f>\, (2.7) and agreement in the graphs of the metric, 



our analysis does seem to indicate that the oscillations observed for this data in fact correspond to the first harmonic 
quasinormal mode of a boson star, however the analysis of the matter field needs further attention. 

Finally, we must remark that we have been unable, using the fundamental and first harmonic modes of an unstable 
boson star, to construct a solution possessing a halo similar to that shown in Figure ||. We do not expect higher 
modes to be of any use here, because the halo is observed to oscillate with the same (single) frequency as the rest of 
the star. Since, as we described at the end of section III, the halo seems to be radiated away over time, we might not 
expect it to be described by the quasinormal modes (which conserve particle number) we have constructed. 
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FIG. 12. First harmonic of an unstable boson star, (a) The solid line shows (j>o5ip\ from the perturbation theory calculations. 
To obtain the squares, we took the simulation data and subtracted the Klein-Gordon field at one instant of time from the data 
at another instant. (The data in the simulations had a spatial resolution four times finer than what is shown in the figure.) 
(b) The squares show the difference between mode obtained from simulation and the mode obtained via perturbation theory, 
scaled relative to the maximum value of 8\(j>\. As we describe in the text, the lack of agreement beyond r ~ 6 may be an artifact 
of simplistic data analysis. The next figure shows that the metric quantities, which depend directly on the matter distribution 
(and thus on \(j>\), show a favorable comparison between the simulations and perturbation theory. 
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FIG. 13. First harmonic of an unstable boson star, (a) The solid line shows the perturbation to a as found from perturbative 
calculations. To plot the squares, we took the simulation data and subtracted the metric function a at one instant of time from 
the data at another instant. (The spatial resolution in the simulation was four times finer than what is shown in the figure.) 
(b) The squares show the difference between the simulation data and the results of linear perturbation theory, scaled relative 
to the maximum value of 8a. The close fit between these results indicates that the oscillations observed in the critical solutions 
correspond to stable oscillatory modes in an unstable boson star. 

VI. HALOS 

We have strong evidence that that the critical solutions correspond to unstable boson stars, but the principal point 
of disagreement is existence of a "halo" of massive field which resides in the "tail" of the solution. It is our contention 
that this halo is not part of the true critical solution, but rather, is an artifact of the collision with the massless field. 

In particular, the halo seems to be a remnant of the original (stable) boson star which is not induced to collapse 
with the rest of the star to form the true critical solution. We find that such a halo is observable in nearly all but the 
most massive (least unstable) critical solutions we have considered, and that its size tends to increase as less massive 
(more unstable) solutions are generated. The fact that the halo thus decreases as we approach the turning point only 
makes sense — a stable boson star very close to the turning point needs very little in the way of a perturbation from 
the massless field to be "popped" over to the unstable branch, and the final, unstable configuration, will, of course, 
be very close to the initial state. 

Additionally, we note that in all cases we have exmained, the field comprising the halo oscillates with nearly the 
same (single) frequency as the rest of the solution. This indicates that the halo is not explainable in terms of additional 
higher-frequency modes. 

As one might expect, the properties of the halo are not universal, i.e. they are quite dependent on the type of initial 
data used. In contrast, the critical solution interior to the halo is largely independent of the form of the initial data. 
To demonstrate this, we use two families of initial data, given by a "gaussian" of Family I in Table | and a "kink" of 
Family II. A series of snapshots from one such pair of evolutions is shown in Figure [l4|. We find different amounts of 
mass transferred from the massless to the massive field for the kink and gaussian data, as shown in Figure [l5|, yet the 
central values of the field oscillate about nearly the same value at nearly the same frequency. Both calculations start 
with identical boson stars with |</>(0,0)| = 0.04 x \/4tt. In the critical regimes, this becomes {\(f>(t, 0)|) = 0.130 x y/4n 
for the solution obtained from the gaussian data, and (\4>(t, 0)|) = 0.135 x \/Att for the kink data. As already noted, 
the oscillation periods are also quite similar, differing by about 3%, and the massess interior to the halo are also quite 
comparable. In particular, it seems quite remarkable that the differences in mass interior to the halo for the two 
families are much smaller than the mass transferred from the real field in either case. 
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FIG. 14. Evolution of r 2 dMc /dr for for two different sets of initial data. Both sets contain the same initial boson star, but 
the massless field <f>z for one set is given by a "gaussian" of Family I (solid line) with ro = 30, and A = 8 whereas for the other 
set <f>3 is given by a "kink" of Family II (dashed line) with ro = 35 and A = 3. The variable A is varied (independently for 
each family) as the parameter p to obtain the critical solution. (Note that after t ~ 60, the massless field has completely left 
the domain shown in the figure.) We have multiplied dMc/dr by r 2 to highlight the dynamics of the halo; thus the main body 
of the solution appears to decrease in size as it moves to lower values of r. The kink data produces a larger and much more 
dynamical halo, but interior to the halo, the two critical solutions match closely — and also match the profile of an unstable 
boson star. Thus, the portion of the solution which is "universal" corresponds to an unstable boson star. 
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FIG. 15. Mc vs. time for the two evolutions shown in Figure |l4|. Mass transfer from the real to the complex field occurs 
from t ~ 30 to t ~ 60, i.e. while the supports of the fields overlap. There is more mass transferred using the kink data, and 
yet the mass falls off rapidly. The mass of the kink data acquires a value very close to the mass of the gaussian data, which is 
itself decreasing slowly with time. We see that, beyond t ~ 250, the difference in mass between the two solutions is very small 
compared with the amount of mass transferred from the real field. 
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If we consider the inner edge of the halo to be where d\<j>\/dr — at some finite radius (e.g., r ~ 5 in Figure g), 
and look at the data between r = and the inner edge of the halo, we find good agreement between this data and 
the profile of a boson star. This can be seen in both Figures || and [l6|. 
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FIG. 16. Mass vs. {\(f>(t,0)\) , the time average of the central value of the field for equilibrium configurations of boson 
stars (solid line), initial data (triangles) and critical solutions (open and filled squares). Arrows are given to help match initial 
data with the corresponding critical solution. Points on the solid line to the left of the maximum mass M ma * ~ 0.633Mj, ; /m 
correspond to stable boson stars, whereas those to the right of the maximum correspond to unstable stars. The data is the same 
as that used for Figure ^, with data from one further evolution added at the bottom of the mass range. The open squares show 
the time average of the mass and \4>{t, 0)| of some critical solutions, and the filled squares show the same quantities evaluated 
between r — and the inner edge of the halo, defined to be the point where d\(j>\/dr = for finite r. The mass of the critical 
solution is in general greater than the mass of the initial data, however the mass inside the halo of the critical solution is less 
than the mass of the initial data. 

We suspect that the halo is radiated over time (via scalar radiation, or "gravitational cooling" ^3|) for all critical 
solutions. We find, however, that the time scale for radiation of the halo is comparable to the time scale for dispersal 
or black hole formation for each (nearly) critical solution we consider. Thus, while we see trends which indicate that 
the halo is indeed radiating, we are not able to demonstrate this conclusively for a variety of scenarios. With higher 
numerical precision, one might be able to more finely tune out the unstable mode, allowing more time to observe the 
behavior of the halo before dispersal or black hole formation occur. 



VII. CONCLUSIONS 



We have shown that it is possible to induce gravitational collapse and, in particular, Type I critical phenomena 
in spherically-symmetric boson stars in the ground state, by means of "perturbations" resulting from gravitational 
interaction with an in-going pulse from a massless real scalar field. Through this interaction, energy is transferred from 
the real to the complex field, and complex field is "driven" and "squeezed" to form a critical solution. The massless 
field is not directly involved in the critical behavior observed in the complex massive field; the critical solution itself 
appears to correspond to a boson star, which, at any finite distance from criticality in parameter space, exhibits a 
superposition of stable and unstable modes. 

Specifically, for initial data consisting of a boson star with nearly the maximum possible mass of M max ~ 
0.633Mp;/m, the resulting critical solution oscillates about a state which has all the features of the correspond- 
ing unstable boson star in the ground state, having the same mass as the initial star. This result is reminiscent of 
the study by Brady et at ||] , who found that the Type I critical solutions for a real massive scalar field corresponded 
to the oscillating soliton stars of Seidel and Suen (JT8). For boson stars with a mass somewhat less than M max , e.g., 
0.9Af max or less, however, we find less than complete agreement between the critical solution and an unstable boson 
star of comparable mass. This is evidenced by the presence of an additional spherical shell or "halo" of matter in the 
critical solution, located in what would be the tail of the corresponding boson star. Interior to this halo, we find that 
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the critical solution compares favorably with the profile of an unstable boson star. Additionally, we have shown that 
the halo details depend on the specifics of the perturbing massless field, and we conjecture that, in the infinite time 
limit, the halo would be radiated away. 

In order to extend the comparison between the critical solutions and boson stars, we have verified and applied the 
linear perturbation analysis presented by Gleiser and Watkins p|, extending their work by providing an algorithm 
to obtain modes with nonzero frequency. We have used this algorithm to give quantitative distributions of mode 
frequency vs. central density of the boson star for the first two modes, as well as to solve for the modes to compare 
with our simulation results. We have found that the unstable mode in the critical solutions have the same growth 
rate as the unstable mode of boson stars, and that the mode shapes also compare quite favorably. We noted that the 
unstable mode of these boson stars was determined much more easily by solving the full nonlinear set of evolution 
equations, rather than via linear perturbation theory. The oscillations observed in the critical solution also indicated 
agreement with first harmonic mode obtained via perturbation theory, however the oscillatory mode in \(/)\ showed 
poor agreement at large radii, and awaits more careful analysis. 

Future work may include simulations of the critical solutions of low mass using higher numerical precision to further 
tune away the initial amplitude of the unstable mode, thus allowing more time to observe the the small halo (i.e., 
whether it is in fact being radiated away). We would also hope to obtain better agreement between simulation 
and perturbation theory for the first harmonic mode of the field \<j>\, perhaps using a more sophisticated method of 
extracting modes from the simulation. Another direction worthy of note would be to begin the simulation with a pulse 
of the complex field (instead of specifically a boson star) tune the height of the pulse to find the critical solutions via 
interpolation, and then compare the resulting critical solutions with our results obtained by perturbing boson stars. 

Finally, we find it worthwhile to investigate similar scenarios for neutron stars. While there have been studies 
regarding the explosion of neutron stars near the minimum mass (e.g., |33j, p4j), we would like to see whether 
neutron stars of non-minimal mass can be driven to explode via dispersal from a critical solution. This may take the 
form of a neutron star approaching the onset of instability via slow accretion, or by being driven across the stability 
graph via violent heating from some other matter source, in a manner similar to the perturbations of boson stars we 
have considered in this paper. 
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APPENDIX A: BOSON STAR MODE FREQUENCIES 



In this Appendix we have tabulated some sample values from the perturbati on th eory c alcula tions. The values and 
uncertainties expressed in the table captions were determined by integrating (4.10) and ( |4 . 1 1 ) to various maximum 
radii, for a range of error tolerances in the integration routines. The values and uncertainties given in the tables were 
chosen to express the variation in our results. 
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TABLE II. 


Shooting Parameters: Fundamental Mode. The values of ( 


/>o(0) are exact. Other quantities 


are given within an 


uncertainty of ±1 in the last significant digit. 
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-1.71E-03 


4.0E-01 
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1.6215E+00 


3.45E+00 


-7.09E-02 


6.0E-01 


1 . 8777E+00 


5.79E+00 


-2.11E-01 



TABLE III. Shooting Parameters: First Harmonic Mode. The values of <^>o(0) are exact, uj is given within an uncertainty of 
±1 in the last significant digit, and the other quantities are given within an uncertainty of ±2 in the last significant digit. 
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APPENDIX B: FINITE DIFFERENCE ALGORITHM 



We approximate the continuum field quantities {a, a, Uj, II2, 1T3, $1, $2, $3, 4>i, §2, ^3} by a set of grid functions, 
quan ti ties w hi ch arc obtained via the solution of finite difference approximations to the partial differential equations 
(2.8), (2.11) - ( |2.14 ) on a domain which has been discretized into a regular mesh (i.e. lattice) with mesh spacing Ar 



in space and At in time. For a grid function u, we denote the value of the grid function in the mesh location j in 
space and n in time by u . 



3 



a. ~ a (nAt, (j ~ l)Ar) , 

where a (nAt, (j — l)Ar) is the corresponding value for the continuum solution. 

The initial data is obtained via "shooting," a standard method of solving ordinary differential equations, in a way 
essentially the same as that found in fl^| . The numerical method used for evolving the system of equations is a 
leapfrog scheme, which is an explicit scheme requiring data at two previous time steps, n and n — 1, to compute a 
value at the next time step n + 1. Given a discretization of scale of order h in time and space, the leapfrog scheme is 
0(h 2 ) accurate. Throughout the mesh, the ratio Acfl = At/ Ar is kept at a constant value, which must be less than 
unity due to the stability requirements of the leapfrog scheme. 

To aid in the presentation of the difference equations, we define the following operators p6[: 
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n+1 



t n 

A u 

3 



2At 



Au. 

3 



3+1 



2Ar 



u = 
- J 



3+1 



Ar 



n n 
r " - 3+1 J-l 



A„u. = 3- 



3 j (r ) 3 -(r ) 3 ' 

V 1 + 1 ; V 1-1 ' 



3 + 

We also define the averaging operator 

r n 1 / n 



^+ j 2 



(n n\ 
U + U. ) , 
J + l 3 J 

which takes precedence over other algebraic operations, e.g. 

r n ( r n\ 1 



h ah 

7 P + 3 



The evolution equations, which are applied to each field Ilj, i — 1, 2, 3} can then be written as: 

n 

A n^= A n(T n ) (Bl) 



03 o V a 



3 



t n .rf r 2 a 



A Q n. =A 3 I — $) -2(000), (B2) 



3 

where the last term in the evolution equation for II is not applied to the massless field. 
Our boundary conditions are as follows: First, by regularity at the origin, we have 

n 

$ x =0 

n+1 

for all n. To obtain II i we employ a "quadratic fit" at the advanced time, 



4n n + l _ n n + l 

n;"~= 2 - 3 , (B3) 



-1 



which is based on the regularity condition, lim r - »o Il(t, r) = Ilo(i) + r 2 Il2(<) + ■ ■ •. 

A significant challenge in the numerical solution of these equations is the problem of the outer boundary condition 
for the massive field. Numerous authors have proposed methods to handle this. Having tried various methods 



including first order expansions of the dispersion relation |18|, sponge filters |27|| , and operator splitting f2q| , we were 



unable to obtain a scheme which produced results superior to the simple Sommerfeld condition one uses for massless 
fields ]29[| . Since, however, the Sommerfeld condition is still inadequate for massive fields, we have chosen to run our 
simulations on a grid large enough that the outer boundary is out of causal contact with the region of interest for 
the time the simulation runs. So, for example, if we are interested in a region < r < 50 and times < t < 400, 
then we place the outer boundary rj > 450. (While unbounded phase velocities are a feature of the Klein-Gordon 
equation, we can argue on physical grounds as well as see quite clearly in simulations that it is the group velocity 
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which is the important quantity in the numerical evolutions, and this is sub-luminal.) Recent work using a shifted 
coordinate system, with a shift vector that is vanishing in some region near r — but increases to unity as r — > rj, 
shows promise as a means of handling the challenge of the boundary condition for the massive field ]3(J], and this 
method may be employed in future work. Thus the outer boundary condition we employ is [j3l| : 

A7 ' A7 ' "J 1 a7 + At 1 (L>4) 

and an analagous equation is used for each II variable. 

After these evolved variables are obtained at the n+l time step, we apply a form of numerical dissipation advocated 
by Kreiss and Oliger p5| . This is applied to both $ and IT in the same manner. So, for instance we set 

n+l n+l 6 / n— 1 n— 1 n— 1 n — 1 n— 1\ 

$ := $ ( $ - 4$ +6$ - 4$ + $ ) , (B5) 

3 3 16 V 3+2 3+1 3 3-1 3-2 7 ' V ; 

where e (0 < e < 1) is an adjustable parameter: typically, we use e = 0.5. 

The preceeding equations describe the "evolution" aspect of the code. The other variables are evolved in a "con- 

n+l n+l 

strained" manner, i.e. they are obtained on the spacelike hypersurface n + l after the fields and II. have 

n+l 

been calculated. The field values c/>, are obtained by updating the value at the outer boundary j — J according to 



. t n 

A 6 

Y J 



n 



and then integrating inward from j = J to j ' = 1 along the spatial hypersurface at n + 1 : 



A.o /''!-,. (B7) 



The Hamiltonian constraint (2.11) can be solved at each time step once all the field variables have been computed 
for the advanced time step. We use the variable A = In a to avoid loss of precision near the origin in the following 
finite difference approximation, which is evaluated at the advanced time step n+l: 



1 - e A r 



A' + A. = M ; + _ [n l + ul+Ill + ^ + ^ + ^ + e A (0? + £)] ) . (B8) 

n+l n+l 

This equation is solved using a pointwise Newton iteration, i.e. given a value of A. (such as A^ — at the 



3 

-1 



origin), we find the next value A. ^ outward along the spatial hypersurface by solving (B8) via Newton's method. 

The slicing condition can be solved once the field variables and the metric function a have been obtained at the 
advanced time step, using the following linear algebraic relation: 



n+l n+l (1/Ar) 

a . . = a . 



3+i 3 (1/Ar)-Z' 



(B9) 



where 

r 

z = n. 



APPENDIX C: DETAILS OF LINEAR STABILITY ANALYSIS 

Following Gleiser and Watkins Jlq] , we write the most general time-dependent, spherically-symmetric metric as 

rfs 2 = - e "frr) dt 2 + e X(t,r) dj ,2 + 

and decompose the complex massive field <p(t, r) via 
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<j>{t,r) = [i>i{t,r) +iip 2 (t,r)]e 



(Cl) 



where tpi and ip2 are real. 

In these variables, the Hamiltonian constraint and slicing condition can be written as 



A' = 



l-e> 



+ r e 



,A-i. 



+ uiP 2 ) 2 + {4>2 - c^i) 2 ] + i>'i + i'i + e A (^i + ^2)) 



(C2) 



1/ = A' + 2- 



2re A K + ^) 



(C3) 



where a prime (') denotes d/dr and an overdot (') denotes d/dt. 
The Klein Gordon equation yields: 



and 



V4' 



- + ^r^) $ + e A (e-^ 2 - l) ^ - e x ~ v $ x + e^^— ^ + ^2) - 2e A -Wa = 
r I I 



2 v' - A' 



'</'2 



e A (e*"u 2 - 1) V2 - e A ">2 + e A ^^-^ 2 - w^i) + 2e A -^i = 0. 



(C4) 



(C5) 



Another equation we will find useful is G\ — 8irGTg , which evaluates to 



2r 2 4 4 



" + + 2^(0^2 - fcVl) + + rt) ) - e- A (V^ + ^) - l$i + 4>1) 



We use equations (|C2|) through (C4) to obtain the equilibrium solutions, by setting 



A(i,r) =A (r) 
i/(t,r) = fo(r) 
V^i(t,r) = 4> (r) 
lM*,r) =0. 



The equilibrium equations are then given by: 



A'o = 



1 - e Ao 



+ r[e A °(a; 2 e-^- 1)^ + ^] 



(C6) 



(C7) 
(C8) 
(C9) 
(CIO) 



(Cll) 
(C12) 



1)^0- 



(C13) 



We now introduce four perturbation fields — 8X(t, r), 8v(t,r), Sipi(t,r) and 8ip2(t,r) — and expand about the equilib- 
rium configuration by writing: 



A(t,r) = X (r) +SX(t,r) 
v(t,r) = vo(r) + 8v(t,r) 

ip2(t,r) = 4> Q (r)8ip 2 {t,r). 



(CI4) 
(CI5) 
(CI6) 
(CI7) 



These last expressions are substituted into (|C2|), (|C3|), (|C4j) and (C6) to obtain the following equations for the 
perturbed quantities: 



23 



(re- Xo 5X)' = r 2 [2(/%6ipi - e~ u ° u 2 (fftSv + 2e~ Va u 2 <t%H\ 
-2e- Vo u4>l8^ 2 + 2e~ Xo (f>' Q ( ( p' S^ 1 + foS^) 



-Ad ^'2 



<Pof>\ 



(C18) 



5v' - 5X' = (u'o -X' + ^jSX- 4re Ao </> 2 1 



Stpi 



(C19) 



+ e A ° (u 2 e~ v ° -l)SX- e x °- Ua uj 2 5v - e^'^S^ - 2e A °" I/ °w<5^ 2 = 



(C20) 



\ 2r 



v'6i>' 



4 

1 



12 - I V 



4 u 4 u y 2 

e- Xo (2<t>' 2 5^ 1 + 2M H'i) + 2^#i] ■ (C21) 
The four equations above can be manipulated such that two variables, bv and Sip2 are eli minat ed, lea ving u s with 



jnly two equations in two unknowns. To obtain the first of these two equations, we subtract (C18) from ( |C20 ) to get 

^-A'\ , 6\' 



Sijji 



2 r 



1-rX'r 



r^ 2 



, e ^o-^ OJ 2_ e \ n 



2e A ° 



1 + e- y °uj 2 



e- X °(^) + rj> cj>> 



(C22) 



To obtai n the other equati on, w e differentiate ( C19| ) with respect to r, and substitute the resulting expression, along 
with ( |C18|) and (|Cl|), into ( |C2l| ) to get 



SX"= --(u' -X' )SX' + 



ajj2 _i_ \" i 2 (^o ~ A o) 2 2z/ o + 
40 o + A + -2 - 



6X 



_Ao — 



dA 



-4 



' to 2 



(C23) 



where, differentiating (Cll) with respect to r we have 



+ r [-v' Q u) 2 e x °-' - ; " 



tl + e^{u 2 e-^ + 1) (X' <g + 2Mo) + 2^ 



(C24) 



(Note that (C22) omits a factor of exp(Ao) which one finds in the ~ 5A/(r 2 (/)g) term of equation (34) in p6|.) For the 
stability analysis, we assume a harmonic time dependence, i.e., 

5itx{t,r) = 5i} 1 {r)e iat 
SX(t,r) = 5X{r)e lat . 



Note that (|C2^ ) and (|C2 4h contain only second derivatives with respect to time. There are good arguments for 
reaT|[|. r- 



assuming a 1 is purely 



161] , so we can determine instability by simply looking for instances where a < 0. 
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As a further consideration, we note that the boson star system admits a conserved Noether current, 



an 



(C25) 



for which the corresponding charge or "particle number" is 



N 



(C26) 



Conventional stability analysis (see, e.g., p2|) demands that we consider only perturbations for which the total 
charge is conserved. Thus we compute the variation in the charge, SN, and work to ensure SN = 0. In practice, 
since we cut off the grid at finite radius, it makes sense to consider the function SN(r), the total charge enclosed in a 
sphere with surface area 4irr 2 . This quantity is 

5N(r) = 




(C27) 



oo. 



where primes denote d/df. (Note that (C27) contains a term involving Sip[, which was not included in equation (35) 
of @.) We then demand that SN -> as r 
The boundary conditions are as follows: 
At r = 0: 



Ao 


= 


^0 


= 


0o 


= 






0o 






1 

~ 3 


SX 


= 


SX' 


= 0. 



3<5A" 



2<j> 2 



1)00 

+ (2{oj 2 + 1) - a 2 ) Sipi 



(C28) 
(C29) 



As 



SN -► 

Siii -> o 
sx -> o. 



To solve the system (C23) and ( 324 ) subject to the above boundary conditions, for a given value of 4>o(0), we resort 
to the method of "shooting," first for the equilibrium solutions, then for the perturbed quantities. Specifically, we 
choose a value for u> and solve the equilibrium equations numerically by integrating outward from r = 0. We do this 
repeatedly, performing a "binary search" on ui (as described in Jl2|) until the boundary conditions for the equilibrium 
quantities are satisfied. 

Due to the linearity of the problem, we can choose Si/ii (0) arbitrarily. We then have two parameters left, namely 
<t 2 and SX"(0). To make matters easy at first, we consider perturbations very close to the transition between stability 
and instability. At the transition point, a 2 is zero. Thus for boson stars near the transition point, we choose a 2 = 
and shoot on the parameter <5A"(0) until the boundary conditions are satisfied. As Gleiser and Watkins jl6| note, the 
transition point occurs at the maximum boson star mass; so we can take two slightly different equilibrium solutions 
near the maximum mass and subtract them to generate solutions which should agree with those obtained from the 
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perturbation problem. We use this method to obtain a trial value of 6X"(0), and also as a way of checking the final 
solution we obtain from the perturbation analysis. 

For more general configurations (a 2 ^ 0), we choose a value of a 2 and shoot on <5A"(0) until we find 6N at the 
outer boundary of the grid to be less than some tolerance value. Then we use the fact (gleaned from experience) that 
if a 2 is too large (too positive), 8N will have a local minimum, the value of which will be less than zero (i.e., 8N(r) 
will dip below zero and then turn back up at larger radii). If a 2 is too low there will be no such local minimum. 
We use these two criteria to select the value of a 2 via a binary search. Thus our two-dimensional eigenvalue-finding 
algorithm consists simply of two (nested) binary searches, one in each direction: For each value of a 2 tried, a full 
binary search on the parameter S\"(0) is performed to drive SN(r max ) — > 0. Then the solution of SN(r) is examined 
for the behavior described above, and a new value of a 2 is selected, and so on until both <5A"(0) and a 2 have been 
found to some desired precision. 
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